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We discuss the GZK horizon of protons and present a method to constrain the injection spectrum of 
ultrahigh energy cosmic rays (UHECRs) from supposedly identified extragalactic sources. This method can 
be applied even when only one or two events per source are observed and is based on the analysis of the 
probability for a given source to populate different energy bins, depending on the actual CR injection spectral 
index. In particular, we show that for a typical source density of 4 x 10 -5 Mpc -3 , a data set of 100 events 
above 6 x 10 19 eV allows one in 97% of all cases to distinguish a source spectrum dN/dE oc E -1 from one 
with E~ 2 ' 7 at 95% confidence level. 
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Introduction — One of the main obstacles to fast 
progress in cosmic ray (CR) physics has been the im- 
possibility to identify individual sources. However, 
there are two pieces of evidence indicating that we are 
at the dawn of "charged particle astronomy." First, 
anisotropies on medium scales have been found combin- 
ing all available data of "old" CR experiments pQ as 
well as in the data from the Pierre Auger Observatory 
(Auger) [2]. Second, the Auger data hint for a correla- 
tion of UHECRs and active galactic nuclei (AGN) [3], 
although this correlations has been contested [4]. Thus 
one may anticipate that the influence of extragalactic 
magnetic fields is small so that UHECRs are not sig- 
nificantly deflected from their initial direction. This 
should be particularly true above the GZK cutoff [5] at 
« 5 x 10 19 eV, when the range of UHECRs is significantly 
reduced by their interactions with photons from the cos- 
mic microwave background (CMB). For instance, for 
typical energy spectra and sources distributed roughly 
homogeneously throughout the universe, 70% of the pro- 
tons with an observed energy of 80 EeV come from 
sources closer than 100 Mpc, even accounting for a 20% 
error in the energy determination. Over such distances, 
the angular spread caused by random magnetic fields 
of 1 nG is typically < 3° for such high-energy protons. 
Deflections in the Galactic magnetic field are expected 
to be of the same order of magnitude [BJ. 

The main reason why no sources have been identi- 
fied yet would be in this scenario that the accumulated 
sky exposure is not yet large enough. While larger ex- 
posures will inevitably increase the number of UHE- 
CRs detected per source, it may take many years until 



enough events are accumulated from even the most in- 
tense source in the sky to allow one drawing a decent in- 
dividual spectrum. The diffuse energy spectrum of CRs 
below E < 4 x 10 19 eV is known with reasonable accu- 
racy and requires a generation spectrum dN/dE oc E~ a 
with a w 2.7 for identical sources — or an appropriate 
distribution of maximal energies E max [7] — while both 
the source and the diffuse spectra at higher energies are 
essentially unknown. It is therefore timely, in the inter- 
mediate phase when sources may be identified by corre- 
lation studies but typically only one or two events per 
source are detected, to ask how the injection spectrum 
can be determined best. 

While first-order Fermi shock acceleration typically 
results in a around 2.1 [H], there exist various models 
that predict either much harder or softer spectra. An 
example for a model with a ~ 1 up to 10 20 eV is the 
acceleration in the electric field around supermassive 
black holes suggested in Ref. [9] [10] that explains also 
the observed properties of large scale jets in AGN [TT] . 
Another possibility to obtain a ~ 1 is to take into ac- 
count a large photon background in the acceleration re- 
gion in the usual shock acceleration [T3]- On the other 
hand, pinch acceleration may serve as an example for 
a = 2.7 [13]. 

In this work, we present an alternative method to set 
constraints on the UHECR source spectrum, suitable 
for the near future of proton astronomy. The basic idea 
to constrain the spectral index of individual sources is 
that, even though the relative weight of different sources 
cannot be known in advance (i.e. before measuring their 
spectra individually), the relative weight of different en- 
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Fig. 1: Distance R in Mpc from which 90% of UHECRs 
arrive with energy > E as function of the threshold en- 
ergy E for £ max = 10 21 eV and a = 2.7. The thin solid 
red line uses CEL in a static Universe as [19], the green 
line uses CEL in an expanding Universe. The blue line 
labeled "SOPHIA" has to be compared to [20]. The red 
line takes into account additionally an experimental en- 
ergy resolution AE/E = 20%. 

ergy bins for a given source is a direct consequence of 
the source spectrum. Now suppose that a minimal en- 
ergy i^min can be identified, above which we can trust 
that the observed CRs come roughly in straight lines 
from their source and, most importantly, sources inside 
the horizon appear with a small enough angular spread 
on the sky that they do not overlap. The energy dis- 
tribution of CRs seen above -E m i n from a given source 
should then reflect the source spectrum (modified by 
the usual propagation effects) , and even if one observes 
only one of them, its energy contains some information 
about the source spectrum. We show how this simple 
argument can be implemented quantitatively for a given 
data set, taking into account UHECR energy losses from 
pure proton sources with supposedly identified distances 
and identical maximum energy. We use this toy model 
to illustrate the basic features of the method and to ex- 
plore its potential power, leaving necessary refinements 
for future work. 

Propagation and horizon scale of UHE protons — In 
Fig. [TJ we show the "90% horizon" - i.e. the distance 
i?90 from which 90% of the UHECRs observed above a 
given energy, E, originate — as function of energy. We 
assume a uniform source distribution with a density 
n s = 4x 10- 5 /Mpc 3 (cf. e.g. Refs. [Hill]) and a power- 
law source spectrum dN/dE oc E~ a with a = 2.7 up 
to the maximal energy E max = 10 21 eV. We used for 
the calculation of photo-pion production the program 
SOPHIA [16 , either taking into account the stochastic- 
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Fig. 2: The distance R in Mpc for which a certain frac- 
tion / of UHECRs arrives with energy > E as function 
of the energy threshold E for 7 = 2.7. From top to bot- 
tom, / = 90% as red line, / = 70% as pink, / = 50% as 
magenta, / = 30% as blue and / = 10% as green line. 



ity of the corresponding energy losses (dotted, blue line) 
or applying the continuous energy loss (CEL) approx- 
imation to its results (dashed, green line). The e + e~ 
pair production losses were taken from Ref. [17] . 

The / = 90% horizon computed within the CEL ap- 
proximation underestimates considerably the full Monte 
Carlo result. The difference increases for a larger "hori- 
zon fraction", / — > 1, and as function of energy for 
E — » -E m ax- There are two reasons for the latter discrep- 
ancy. First, the energy transfer per interaction, y, in- 
creases with energy and violates more and more strongly 
the formal requirement y <C 1 needed for the applicabil- 
ity of the CEL approximation. Second, the flux taking 
into account the stochastic nature of the energy losses in 
pion production remains finite for E — > E max , while in 
the CEL approximation no particles with E = E max can 
reach the observer from a source at a finite distance [TS] • 

In a realistic experiment, the primary energy can 
only be reconstructed with a finite precision. Assum- 
ing a Gaussian (in log E) experimental uncertainty of 
AE/E = 20%, we computed the 90% horizon as a func- 
tion of the measured CR energy, for the same conditions 
as above. The two resulting curves are also shown in 
Fig. Q] Since the CR spectrum is falling steeply, the 
misinterpretation of lower energy events as high energy 
ones has a larger impact than the reverse, which in turn 
leads to an increase of the estimated horizon scale. At 
low energies, say <5x 10 19 eV, the observed spectrum 
approximates well to a power-law and the energy res- 
olution only affects the absolute flux, not the relative 
fluxes relevant for Rgo(E). 
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The horizon scale for UHE protons and nuclei was 
recently discussed also in Refs. [TH] [2D]- In Fig- Q] we 
compare our calculations to those of Refs. [19] and [20] 
for proton primaries. In Ref. [19j . Harari et al. pre- 
sented results (shown as orange line) using the CEL ap- 
proximation and assuming a static Universe, our result 
for the same assumptions is shown with a thin solid red 
line. Both calculations agree well at moderate energies 
E = 80 — 100 EeV, while there is some disagreement 
both at high and low energies. However, the differences 
at low energies between the two calculations are much 
smaller than the differences between those calculations 
and the more correct CEL calculation in the ACDM 
model for the expanding Universe, presented with a 
green line. 

All results using the CEL approximation differ in 
shape as a function of energy from the calculationss us- 
ing SOPHIA for pion production either directly (blue 
line), or using the SOPHIA results in a kinetic equation 
approach as in Ref. [20] (magenta line). The agreement 
between the latter two results is almost perfect at all 
energies. 

As an illustration, we show in Fig. [2] the horizon dis- 
tance corresponding to different CR fractions. Specifi- 
cally, we plot the distance Rf below which a given frac- 
tion / of the UHECRs reach the Earth with an energy 
larger than E, as a function of that energy, for / = 10%, 
30%, 50%, 70% and 90% (using always SOPHIA and 
AE/E = 20%). 

Estimation of the spectral index — Since the angu- 
lar resolution of cosmic ray experiments is poor by as- 
tronomical standards, the identification of individual 
sources requires a relatively large angular distance be- 
tween them. This can only hold for sufficiently high 
energies such that the horizon scale is small, say of the 
order of 100 Mpc, leaving a limited number of sources 
over the sky. Defining as horizon, within which 90% of 
all CRs observed above a given energy were emitted, we 
find from Fig. [5] that a horizon of 100 Mpc corresponds 
to a threshold energy of E = 1 X 10 20 eV. At present, 
the importance of deflections in extragalactic magnetic 
fields above this energy is unclear. As soon as sources 
are detected, one will be able to set an upper limit and to 
a certain extent reconstruct the extragalactic magnetic 
field. Here, we limit ourselves to the optimistic scenario 
where deflections in extragalactic magnetic fields are not 
much larger than the combined effects of the Galactic 
magnetic field and the experimental angular resolution. 

At present, the picture of uniformly distributed, ex- 
tragalactic UHECR sources having all the same lumi- 
nosity and the same injection spectrum is able to de- 
scribe well the observed energy spectrum in a broad 



energy range from a fewxlO 17 eV or a fewxlO 18 eV up 
to the GZK cutoff, depending on the assumed source 
composition [2T1 122j . 

We first produce a Monte Carlo (MC) sample by gen- 
erating sources with constant comoving density n s = 
4 x 10~ 5 Mpc -3 up to a maximal redshift of z = 0.1. 
Then we choose a source i according to the declina- 
tion dependent exposure of Auger, with an additional 
weight chosen according to the source distance. Fi- 
nally, we generate a CR with an initial energy drawn 
randomly according to the assumed injection spectrum, 
dN/dE oc E~ a " , and propagate it until it either reaches 
the Earth distance or loses energy down to below E m [ n . 
In the former case, we then apply an energy-dependent 
angular deflection to mimic the effect of the Galactic 
magnetic field, with a shift perpendicular to the Galac- 
tic plane equal to 5b = 2°(£:/10 20 eV)~ 1 , where this 
magnitude is motivated by the results of Ref. [6] . The 
chosen magnetic field likely overestimates deflections far- 
away from the galactic plane in most of models. How- 
ever, we consider this choice as a conservative upper 
limit. Finally, we deflect the CR direction to account 
for a finite experimental angular resolution, taking the 
Auger surface detector as a reference [23] , with a spher- 
ical Gaussian density oc exp(— £ 2 /(2af)) sm(£)d£, where 
ai — 0.85° and I is the angular distance. 

After having generated N cosmic rays, we perform 
a correlation analysis between the CRs and the sources. 
First, we identify as "the source" of a given CR the 
source with the smallest angular distance I to the ob- 
served CR arrival direction and maximal distance R — 
100 Mpc. Inside this region, there are around ~ 160 
sources for chosen density n s — 4 x 10~ 5 Mpc~ 3 . Such 
a small number makes the probability negligible that 
sources overlap, if they are uniformly distributed. This 
probability increases, if sources follow — as expected — 
the large-scale structure of matter and may constitute a 
real limitation to resolve single sources in cluster cores. 

Additionally, we require that the angular distance I 
be smaller than a prescribed value, £ max - Next, having 
pre-defined an energy E 2 that divides the whole energy 
range into two large bins, we count for each source i the 
numbers N^i and AT ii2 of high energy (E > E 2 ) and low 
energy events (E m i n < E < E 2 ), respectively. Given the 
corresponding fractions fi(a) and /2(a) = 1 — fi{a) of 
Ni = Ni,i + Ni : 2 events expected from a source at the 
identified distance for an arbitrary value of the spectral 
index a, we calculate with a binomial distribution the 
probability, 

Pi{N i>1 ,N ii2 \a)= ' , ', A (")/ 2 (<*)> (1) 

J M 1 liVi O' 



4 



that the observed numbers N^j are consistent with the 
value ao used in the MC. Considered as a function of 
a, this probability distribution has the true value ao as 
its expectation value, if our procedure is unbiased, and 
measures how strongly the data disfavor a differently 
assumed value q^qo. 

Since the different sources emit CRs independently 
from one another, we can simply multiply the single 
source probabilities Pi(N iy i, Ni^\a) to obtain the global 
probability of a given data set with N s identified sources: 

KW,i}i=i,;v» = Y[pi(fi,i,fi,2\a)- (2) 
»=i 

The basic outcome of a sample of MC simulations 
for fixed parameters 9 — a>o . . . is thus a binned distribu- 
tion, f(p\0), giving the fraction / of MCs producing the 
value p. With how much confidence can we distinguish 
these distributions for two different B\ and 6*2? Clearly, 
the smaller the overlap of the two distributions, the eas- 
ier the two parameter sets Qi can be distinguished. 
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Fig. 3: Distribution of probability of reconstructed 
power law spectrum if real power law spectrum is 
a — 2.7, angle £ max = 4°. In all cases -Emm = 60 EeV. 
The red line is for ot — 1.1 and the blue line for a = 2.7. 

We study now the possibility to distinguish differ- 
ent values of the injection spectrum of CRs in more 
detail. As simplifying assumption we assume that the 
injection spectrum of all sources is the same, i.e. in par- 
ticular that the maximal energy of all sources is iden- 
tical. This assumption allows us to study the spectra 
only above ~ 4 x 10 19 eV, because at lower energies a 
spectral index a < 2.6 requires either additional Galac- 
tic sources or a non-uniform source distribution. In the 
latter case, either the source density or the luminos- 
ity of single sources should increase as function of red- 
shift, n(z) = n (l + z) m and L(z) = L (l + z) m respec- 



tively, or the maximal energy of sources is distributed 
as dn/dE max oc a [TJ. Moreover, we consider only 
two extreme cases, namely a power- law with ao = 1.1 
and ao = 2.7. 

In Fig. [3] we compare the distributions of probabili- 
ties obtained from Eq. choosing as true value a — 
2.7, as source density as always n s = 4 x 10~ 5 /Mpc 3 , 
as number of CRs N = 100, and £ max = 4°. The red 
solid line is the distribution of probabilities obtained as- 
suming a — 1.1, while the blue dashed line corresponds 
a = 2.7. The two curves have only a small overlap, since 
the probabilities using the correct a are rather narrowly 
concentrated around p = 1 , while the probability distri- 
bution using the wrong a extends from extremely low 
values up to one. Thus an experimental differentiation 
between different injection spectra seems possible, even 
if only one or, in few cases, two events per source are 
detected, as it is the case for the chosen parameters in 
Fig. [3] This constitutes the main result of our work. 

We quantify the chances to distinguish two differ- 
ent spectral indices in the following way: We calculate 
the area A corresponding to the desired confidence level 
(C.L.), A, starting from 1 to the left using the best- 
fit distribution (e.g. the blue line in Fig. [3j and obtain 
thereby as its lower boundary pa- Thus only in 1 — pa 
cases we will obtain by chance a lower probability using 
the correct test hypothesis. Next we count how large is 
the area B of the wrong test hypothesis on the left of 
Pa- As final answer we obtain that in the fraction B of 
all cases we can distinguish between the two hypotheses 
with C.L. A. 

Let us illustrate this procedure for the case consid- 
ered above, choosing as confidence level A = 95%. The 
green dashed-dotted vertical line in Fig. [3] enclosing 95% 
of the area of the true (blue) distribution determines 
P95 = 0.056. The area of the red curve on the left of 
Pq5 — 0.056 is B = 0.971. Hence one can exclude in 
B = 97.1% of cases with at least 95% C.L. the expo- 
nent a = 1.1 for the spectrum, if the true exponent is 
a a = 2.7. 

In addition to the rather extreme cases of the spec- 
tral indices above, we investigated the ability of the 
method to distinguish between any of them and an in- 
termediate value of «o = 2.0, often considered in the 
context of astrophysical particle acceleration. As an il- 
lustration, we found that with a data set of 100 cosmic 
rays above 6 x 10 19 eV, it is possible in 50% of the cases 
to discriminate ao = 2.0 from a value of either 1.1 or 2.7 
with a C.L. of 95%. Likewise, for a data set of 200 cos- 
mic rays above 4 x 10 19 eV (i.e. for essentially the same 
exposure of the sky, but with a lower energy threshold) , 
an injection spectral index of 1.1 can be discriminated 
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against ao = 2.0 with a C.L. of 95% in 90% of the cases, 
while an injection spectral index of 2.7 can be discrim- 
inated against cto = 2.0 with a C.L. of 95% in 70% of 
the cases. 

Summary — We have proposed a method to estimate 
the generation spectrum of individual extragalactic CR 
sources that is well-suited for the time when only one 
or two events per source are detected. An important in- 
gredient of this method is the relative fraction of events 
contained in a prescribed energy interval. Therefore we 
have recalculated the horizon scale of ultra-high energy 
protons, taking into account a reasonable energy reso- 
lution, similar to that of Auger. 

We have demonstrated for a toy-model the potential 
of this method, finding that around 100 events above 
6 x 10 19 eV are required to distinguish with 97% prob- 
ability at least at the 95% C.L. the two extreme cases 
a = 1.1 and 2.7. A differentiation between a's that are 
more similar will be clearly more challenging. An injec- 
tion spectral index of 2.0 can still be distinguished from 
the two above values with a 95% C.L. in the majority 
of cases (with the same statistics). 

Several of the issues we have neglected, like the effect 
of a possible -E max distribution, should be included in a 
more complete study as soon as experimental data will 
be available. A proper estimation of a also requires to 
quantify the bias introduced e.g. by misidentified events. 
In general, it proves more efficient to remove from the 
data set the doubtful events (e.g. in regions where a 
given catalogue used to identify sources is known to 
be incomplete, or when several sources at different dis- 
tances are identified over a small region of the sky, with 
possible overlap due to magnetic deflection or poor an- 
gular resolution), and apply the method with a corre- 
spondingly smaller statistics. Sources physically clus- 
tered in the universe are not a problem here, since they 
are located essentially at the same distance from the 
Earth and thus suffer from the same attenuation during 
propagation. 
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